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Предложен способ аппроксимации плотностей статистических распределений полиномами Чебышева-Эрмита для дальнейшей 
генерации последовательностей случайных чисел, закон распределения которых близок к аппроксимируемому закону. Разрабо- 
тана методика регуляризации решения системы линейных алгебраических уравнений, возникающей при решении задачи ап- 
проксимации плотности распределения. Исследованы возможности применения предложенного способа для аппроксимации 
гладких и островершинных распределений, имеющих одну моду, а также двухмодальных распределений. 


1. Введение 

При моделировании поведения какого-либо 
процесса или объекта, находящегося под воздей- 
ствием случайных факторов, возникает задача по- 
лучения случайных чисел, закон распределения ко- 
торых близок к закону распределения, наблюдае- 
мому в эксперименте [1]. В этом случае зачастую не 
важно знать тип и значения параметров закона ра- 
спределения генерируемых случайных чисел, важ- 
но лишь, чтобы этот закон соответствовал закону 
распределения, полученному при наблюдении за 
поведением исследуемого объекта или процесса. 
Один из подходов к решению указанной задачи за- 
ключается в нахождении некоторой аппроксими- 
рующей функции, которая описывает эксперимен- 
тальный закон распределения с требуемой точно- 
стью. На основе этой функции искомые случайные 
числа могут быть получены из равномерных слу- 
чайных чисел методом обратной функции [2]. 

Несмотря на кажущуюся простоту такого под- 
хода, при его практической реализации возникает 
ряд проблем, связанных, в первую очередь, с тре- 
бованиями, которые предъявляются к функциям, 
аппроксимирующим плотность распределения и 
функцию распределения. К этим требованиям от- 
носятся: 1) положительность функции, аппрокси- 
мирующей плотность распределения, на всей ее 
области определения и 2) равенство интеграла, взя- 


того от аппроксимирующей функции по всей обла- 
сти ее определения, единице. Кроме того, зачастую 
возникают проблемы, связанные с низкой устой- 
чивостью получаемых решений. Одним из спосо- 
бов, позволяющих успешно решить задачу аппрок- 
симации плотности распределения, является под- 
ход, основанный на применении ортогональных 
полиномов Чебышева-Эрмита [3]. 


2. Постановка и решение задачи аппроксимации 
плотности распределения полиномами 
Чебышева-Эрмита 

Как известно, многочлены Чебышева-Эрмита 
Н(х) ортогональны на всей числовой оси х с весо- 
вой функцией 

Н(х) = е~ х \ (1) 


и, соответственно, удовлетворяют соотношениям 


| Н і {х)Н ] (х)Н(х)сіх 
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В случае, когда используются ортонормирован- 
ные многочлены Чебышева-Эрмита, значение пе- 
ременной ^ равно единице. 

Вычисление многочленов Чебышева-Эрмита 
может осуществляться при помощи формулы Ро- 
дрига, имеющей вид [3]: 
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я,Тх) = (-іуѵѴ 2 ) < " ) , 

либо при помощи рекуррентного соотношения вида 
Н п+1 (х) = 2хН п (х)-2пН п _ 1 (х). 

Как показано в работе [3], частичные суммы ря- 
дов Фурье для некоторой функции у(х) по ортого- 
нальным многочленам являются наилучшими при- 
ближениями данной функции на конечном интерва- 
ле [ а ; Ь]. Тот факт, что весовая функция для ортого- 
нальных многочленов Чебышева- Эрмита является 
решением известного уравнения Пирсона [3], делает 
предпочтительным использование данного семей- 
ства многочленов для аппроксимации плотностей ра- 
спределений с учетом того, что они ортогональны на 
всей числовой оси. Следует отметить также, что при 
разложении функции у(х) в ряд Фурье по многочле- 
нам Чебышева- Эрмита к ней предъявляются мини- 
мальные требования в смысле ее интегрируемости. 

В рамках данной работы аппроксимация плот- 
ности распределения, наблюдаемой в эксперимен- 
те, производилась при помощи функции следую- 
щего вида: 

У(х) = X а і Н і (*)> (2) 

2=0 


где п — максимальный порядок многочлена; а { - 
постоянные коэффициенты. 


В этом случае задача аппроксимации сводится к 
нахождению таких значений коэффициентов а ь 
при которых выполняется условие 
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где т - число интервалов гистограммы, построен- 
ной по экспериментальной выборке. 

Дифференцируя выражение (3) по неизвестным 
коэффициентам а ь получаем систему из п линей- 
ных уравнений следующего вида: 
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Н к (х . ), к - 0...п. (4) 


После ряда очевидных преобразований система 
уравнений (4) может быть представлена в следую- 
щей матричной форме: 

Ва = у, (5) 


где а - вектор искомых коэффициентов а ь а эл- 
ементы Ь и матрицы В и компоненты у { вектора у 
вычисляются следующим образом: 

т 

\/=Х я ,ыды, (6) 

к = 1 


у,=Ха я ,ю- о 

к = 1 

Очевидно, что матрица В является симметрич- 
ной, и ее элементы, расположенные на главной ди- 
агонали, отличны от нуля. Следует отметить также, 
что элементы матрицы В, не расположенные на 


главной диагонали, также в общем случае отличны 
от нуля, что обусловлено следующими факторами: 

1. в выражении (6) при суммировании полиномов 
не используется весовая функция (1); 

2. интервал значений переменной х, на котором 
производится суммирование, конечен; 

3. суммирование производится только в конечном 
множестве точек х ь к=1,...,т. 

Решая систему (5), можно найти вектор иско- 
мых коэффициентов а г Предложенный подход по- 
зволяет с необходимой точностью аппроксимиро- 
вать плотности одномодальных распределений, 
форма которых близка к форме нормального ра- 
спределения (рис. 1). 
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Рис. 1. Аппроксимация плотности р(х) нормального распре- 
деления при помощи полиномов Чебышева -Эрмита 

Однако осуществить аппроксимацию плотно- 
стей распределений, форма которых значительно 
отличается от нормального распределения, в част- 
ности двухмодальных распределений, не удается. В 
частности, наблюдается появление «всплесков» ап- 
проксимирующей функции. Кроме того, на краях 
области определения в некоторых случаях возни- 
кают участки, где аппроксимирующая функция 
опускается ниже оси абсцисс. Примером такого 
поведения описанного алгоритма являются резуль- 
таты аппроксимации плотности двухмодального 
распределения (рис. 2). 

Как показано на рис. 2, варьирование макси- 
мального порядка аппроксимирующего полинома 
не позволяет осуществить аппроксимацию таким 
образом, чтобы выполнялись сформулированные 
выше требования. Данный факт обусловлен, по- 
видимому, слабой устойчивостью решения систе- 
мы (5) к ошибкам счета, возникающих при вычи- 
слении полиномов, входящих в сумму (2), причем, 
очевидно, между порядком этих полиномов и вели- 
чинами ошибок имеется некоторая взаимосвязь. 

Для повышения устойчивости решения систе- 
мы уравнений (5) авторами был применен широко 
известный метод регуляризации, предложенный 
А.Н. Тихоновым [4]. В рамках данного метода точ- 
ное решение а системы уравнений (5) заменяется 
некоторым приближенным решением другой си- 
стемы, которая имеет следующий вид: 

(В + у Е)а = у, у >0, (8) 

где Е - единичная матрица, у- некоторый малый ко- 
эффициент, называемый параметром регуляризации. 
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в) л =15 г) п = 20 

Рис. 2. Результаты аппроксимации плотности р(х) двухмодального распределения при различных значениях порядка аппрок- 
симирующего полинома п 


Рассматривая эту систему, можно непосред- 
ственно видеть, что при: 1) больших значениях у 
решение ее сильно отличается от решения системы 
(5), однако является устойчивым; 2) малых /реше- 
ние близко к решению системы (5), однако устой- 
чивость его уменьшается. Следовательно, необхо- 
димо выбирать такое значение у, при котором точ- 
ность и устойчивость получаемого решения явля- 
ются оптимальными [4]. 

Как показали проведенные исследования, в об- 
щем случае невозможно подобрать такое значение 
у, при котором обеспечивалась бы достаточная 
точность аппроксимации при отсутствии колеба- 
ний аппроксимирующей кривой на краях ее обла- 
сти определения. Данное явление связано, очевид- 
но, с существенно различными значениями оши- 
бок счета, возникающих при вычислении полино- 
мов различных степеней, и, как следствие, с раз- 
личными значениями ошибок задания элементов 
матрицы В и вектора у в системе (5). Таким обра- 
зом, представляется логичным заменить в (8) еди- 
ничную матрицу Е на некоторую матрицу весовых 
коэффициентов М, сформированную на основа- 
нии имеющейся информации о величинах ошибок 
вычисления каждого из полиномов, либо хотя бы о 
соотношении значений этих ошибок. При этом си- 
стема (8) преобразуется к следующему виду: 

(В + аМ)а = у, а > 0, (9) 

где а - малый положительный параметр. 

Известно, что ошибки вычислений на ЭВМ, 
проводимых с использованием чисел с плавающей 
точкой, определяются множеством факторов, влия- 
ние которых в большинстве практических случаев 
учесть невозможно [5]. В качестве примеров подоб- 
ных факторов можно привести особенности архи- 
тектуры центрального процессора, используемой 


операционной системы, языка программирования, 
версии компилятора и т. д. В связи с этим одним из 
возможных способов формирования матрицы М, 
входящей в систему уравнений (9), является анализ 
значений элементов матрицы В и вектора у. 

Пусть каждый из полиномов Щх) вычисляется 
в точке х к с ошибкой А Н[. Тогда выражения (6, 7) 
можно переписать следующим образом: 

т 

Кі = X (н , <Л ) + А К X я , (Л ) + Щ ) = 

к = 1 

( 10 ) 


Уі = X А ( н і (к )+Ю- (Н) 

к = 1 

Разделив обе части равенства (10) на соответ- 
ствующие части равенства (11), пренебрегая в вы- 
ражении (10) членом высшего порядка малости, 
получаем: 

X (Я, (х к )Н 1 (х к ) + Я, (х к )М- к + Н і (х к )Ак [ ) 

(к = 1=1 .(12) 

у ‘ ІуЛЛ(Щ+К) 

к = 1 


=х 


Г НЩЦНЦхЦ + Н Цх к )Ыі[ + Л 


к = 1 


у 


V 


+Я / (х*)Д#+ДАда 


У 


При і=і выражение (12) преобразуется к сле- 
дующему виду: 


т 

^(Н і (х к )(Н : (х к ) + 2А/()) 

°І,І _ к = 1 


т 

^уЦнщп+К) 


) 

к= 1 

т 

X -Й ^ (у ) 

к = 1 



Важными особенностями отношения (13) явля- 
ются следующие: 
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в) п = 1 5 г) п = 20 

Рис. 3. Влияние параметра а на результаты аппроксимации плотности р(х) двухмодального распределения при различных 
значениях п 


1 . при стремлении ошибки А Н[ к нулю отношение 
(13) стремится к отношению точных значений 
элементов Ь ц и у ( ; 

2. в зависимости от знака ошибки отношение (13) 
может отличаться от отношения истинных зна- 
чений элементов Ъ ц и у і как в меньшую, так и в 
большую сторону; 

3. примерное равенство в (13) справедливо только 
при і =у, т. к. значения ошибок вычисления по- 
линомов различных степеней могут сильно раз- 
личаться. 


На основании сказанного выше, сформируем 
матрицу М следующим образом: 


м = {щ.і } = < 


Ъ ч 


У; 

0, і Ф і 


. і=] 


Определим теперь значение малого коэффици- 
ента а, входящего в выражение (7). Для удобства 
его определения пронормируем матрицу М следу- 
ющим образом: 


у, увеличение параметра а не столь негативно ска- 
зывается на точности аппроксимации. При а> \ на- 
блюдается подавление колебаний аппроксимирую- 
щей кривой, и вместе с тем точность аппроксима- 
ции снижается незначительно. В целом, наилуч- 
шие результаты аппроксимации в смысле макси- 
мальной точности и минимальных колебаний были 
получены при значениях ае(2;5). 

Важным достоинством предложенного алгорит- 
ма является высокая устойчивость получаемых ре- 
шений при увеличении порядка п полинома (2). На 
рис. 4 приведены аппроксимирующие кривые, по- 
лученные при п =50. 



- 2-10123 
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м = {щ.і } = 


Щ,і 


тах{т і . ) 
0, і Ф у 


> * =У 


Таким образом, значения элементов матрицы М 
будут находиться в интервале от нуля до единицы. 

Рассмотрим влияние значения параметра а на 
результаты аппроксимации при различных поряд- 
ках полинома п (рис. 3). 

Как видно из рис. 3, варьируя значение а, мож- 
но добиться достаточной точности аппроксимации 
при отсутствии «всплесков» и отрицательных зна- 
чений аппроксимирующей кривой. Необходимо 
отметить, что, по сравнению с влиянием параметра 


Рис. 4. Результаты аппроксимации плотности р(х) двухмо- 
дального распределения при различных значениях 
параметра сс и п = 50 

Как показывают результаты, полученные при 
аппроксимации экспериментальной двухмодальной 
плотности распределения, при увеличении п точ- 
ность аппроксимации увеличивается (рис. 5). Без 
использования регуляризации точность аппрокси- 
мации максимальна, однако с ростом п наблюдают- 
ся существенные колебания аппроксимирующей 
функции на краях области ее определения. Приме- 
нение метода регуляризации А.Н. Тихонова [4] по- 
зволяет подавить колебания аппроксимирующей 
функции, но точность аппроксимации при этом 
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снюхается в 5. ..6 раз. Применение изложенного в 
данной работе алгоритма регуляризации позволяет 
подавить колебания аппроксимирующей функции и 
получить точность аппроксимации, сравнимую с 
той, которая имеет место без использования регуля- 
ризации. Таким образом, результаты, приведенные 
на рис. 4 и 5, являются подтверждением высокой 
устойчивости решений, получаемых с помощью 
разработанного алгоритма, к ошибкам вычисления 
полиномов Чебышева- Эрмита высоких степеней. 



п 


Рис. 5. Зависимости среднеквадратического отклонения 5 от 
порядка полинома п, построенные при аппроксима- 
ции плотности распределения, наблюдаемого в экспе- 
рименте, с использованием регуляризации и без нее 


3. Применение результатов аппроксимации 
плотностей распределений для получения 
последовательностей случайных чисел 
с заданным законом распределения 

Замечательной особенностью разработанного 
способа аппроксимации плотностей распределений, 
основанного на применении полиномов Чебышева- 


Эрмита, является возможность аналитического вы- 
числения интеграла от аппроксимирующей функ- 
ции. В работе [3] показано, что определенный инте- 
грал от полинома Чебышева- Эрмита к-то порядка 
Н к (х) может быть представлен следующим образом: 



Н к+1 (Ь)-Н к+1 (а) 
2 (к + 1) 


Следовательно, определенный интеграл от ап- 
проксимирующей функции (2), взятый на интерва- 
ле |Х; х 2 ], может быть вычислен при помощи про- 
стого выражения: 


Х 2 п 

I у(.х)сІх = ^ а , 

т. 2=0 


Н і+ 1 (х 2 ) Н і+1 (х 1 ) 
2(7 + 1) 



Таким образом, для вычисления интеграла от 
аппроксимирующей функции не требуется приме- 
нения каких-либо численных методов. Отметим, 
что, поскольку функция (2) аппроксимирует плот- 
ность некоторого распределения, то интеграл от 
(2), вычисленный при помощи выражения (14), ап- 
проксимирует, очевидно, функцию данного ра- 
спределения. Эта особенность позволяет при гене- 
рации случайных чисел, закон распределения ко- 
торых близок к экспериментальному закону ра- 
спределения, использовать высокопроизводитель- 
ные методы вычисления обратной функции, такие 
как метод дихотомии, метод Ньютона и др. 

На рис. 6 приведены гистограммы, аппрокси- 
мирующие их функции, а также интегралы от ап- 
проксимирующих функций, построенные для мо- 
дельных выборок и выборок значений случайных 
величин, наблюдаемых в эксперименте. 

Как видно из рис. 6, удовлетворительные ре- 
зультаты аппроксимации для различных законов 
распределения получаются при существенно раз- 



а) п = 8; 5= 0,003 



О 0,5 

в) п = 35; 5= 0,010 


X 


б) п = 10; 5 =0,021 



г) п= 33; 5 1 = 0,020 


Рис. 6. Г истограммы, аппроксимирующие их функции, а также интегралы от этих функций; а) нормальное распределение; 
б) гамма-распределение; в) бета-распределение; г) распределение Коши; сс=2 
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личающихся порядках полинома п. Кроме того, ап- 
проксимация островершинных распределений, та- 
ких как распределение Коши, не может быть осу- 
ществлена с достаточно высокой точностью. Уве- 
личение п в подобных случаях не приводит к улуч- 
шению качества аппроксимации, т. к. при этом 
увеличиваются колебания аппроксимирующей 
функции. Результаты, приведенные на рис. 6, по- 
зволяют сделать вывод о том, что область примене- 
ния разработанного подхода ограничена; наилуч- 
шим образом данный подход может быть использо- 
ван при аппроксимации достаточно гладких плот- 
ностей распределений. 

Несмотря на этот факт, интеграл от аппрокси- 
мирующей функции (2) близок к интегралу от ап- 
проксимируемой плотности распределения. Таким 
образом, возможно применение разработанного 
подхода при генерации случайных чисел, подчи- 
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няющихся требуемому закону распределения, ме- 
тодом обратной функции. 

Заключение 

Разработанный способ аппроксимации плотно- 
стей распределений на основе ортогональных по- 
линомов Чебышева- Эрмита может быть с успехом 
применен при обработке экспериментальных вы- 
борок и перспективен при решении широкого кру- 
га научных и практических задач, связанных с мо- 
делированием случайных величин. Предложенная 
методика регуляризации позволила получить 
устойчивые решения даже при достаточно высоких 
порядках (я<50) аппроксимирующих полиномов. 
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